ICCB Environmental data download
In this script we are downloading current and future environmental data (as rasters) to use in species distribution models (SDMs) for koalas (Phascolarctos cinereus) in the South-East Queensland (SEQ) region.
Import packages
Load South East Queensland (SEQ) boundary
We start by defining our study area, which is the South East Queensland (SEQ) region. We will use the Local Government Areas (LGA) shapefile to define the extent of SEQ.
https://qldspatial.information.qld.gov.au/catalogue/custom/detail.page?fid={3F3DBD69-647B-4833-B0A5-CC43D5E70699}
Code
Reading layer `Local_Government_Areas' from data source
`/Users/scottforrest/Library/CloudStorage/OneDrive-QueenslandUniversityofTechnology/PhD - Scott Forrest/GIT/ICCB_geospatial_tools_conservation/Session 3/Data/Environmental_variables/Local_Government_Areas.shp'
using driver `ESRI Shapefile'
Simple feature collection with 78 features and 8 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 137.9946 ymin: -29.17927 xmax: 153.5519 ymax: -9.087991
Geodetic CRS: GDA94
Coordinate Reference System:
User input: GDA94
wkt:
GEOGCRS["GDA94",
DATUM["Geocentric Datum of Australia 1994",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore."],
BBOX[-60.55,93.41,-8.47,173.34]],
ID["EPSG",4283]]
Code
# Convert to WGS84
LGA <- LGA %>% st_transform(7856)
# Select local govt. areas for South East Queensland
LGA_SEQ <- LGA %>%
filter(lga %in% c("Brisbane City",
"Moreton Bay City",
"Logan City",
"Ipswich City",
"Redland City",
"Scenic Rim Regional",
"Somerset Regional",
"Lockyer Valley Regional",
"Gold Coast City",
"Sunshine Coast Regional",
"Toowoomba Regional",
"Noosa Shire"))
ggplot() +
geom_sf(data = LGA, color = "black") +
geom_sf(data = LGA_SEQ, fill = "purple3", alpha = 0.5, color = "black", size = 0.2) +
theme_minimal() +
theme(legend.position = "none") +
labs(title = "Local Government Areas Queensland (SEQ in purple)")Code
Merge into a single polygon
Save SEQ extent for other scripts
Load current environmental data
Layers were made available to us by the EcoCommons team and were created by Toombs and Ma (2025):
Toombs, N., and Ma S., 2025, A High-Resolution Dataset of 19 Bioclimatic Indices over Australia, Climate Projections and Services – Queensland Treasury, Brisbane, Queensland. [https://longpaddock.qld.gov.au/qld-future-climate/data-info/tern/]
Code
files <- list.files("Data/Environmental_variables/Current_climate_QLD",
pattern = ".tif$",
full.names = TRUE)
# Load all bioclim rasters
current_bioclim <- lapply(files, terra::rast)
# Make into one raster stack
current_bioclim <- rast(current_bioclim)
# Plot the current bioclimatic variables
plot(current_bioclim)class : SpatRaster
dimensions : 681, 841, 19 (nrow, ncol, nlyr)
resolution : 0.05, 0.05 (x, y)
extent : 111.975, 154.025, -44.025, -9.975 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
sources : bioclim_01.tif
bioclim_02.tif
bioclim_03.tif
... and 16 more sources
names : bioclim_01, bioclim_02, bioclim_03, bioclim_04, bioclim_05, bioclim_06, ...
[1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n MEMBER[\"World Geodetic System 1984 (G2296)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"
Code
class : SpatRaster
dimensions : 834, 898, 19 (nrow, ncol, nlyr)
resolution : 5590.925, 5590.925 (x, y)
extent : -4407995, 612655.7, 4234346, 8897178 (xmin, xmax, ymin, ymax)
coord. ref. : GDA2020 / MGA zone 56 (EPSG:7856)
source(s) : memory
names : bioclim_01, bioclim_02, bioclim_03, bioclim_04, bioclim_05, bioclim_06, ...
min values : 5.182303, 5.520304, 35.25507, 105.7126, 16.82971, -4.954343, ...
max values : 29.453394, 16.870607, 61.84691, 680.0889, 42.45033, 22.998774, ...
Mask and crop to SEQ extent
Code
Load future environmental data
Here we load outputs from a moderate-high emissions shared socio-economic path scenario (SSP 3.70) for the year 2090 (2080 - 2099).
Code
files <- list.files("Data/Environmental_variables/Future_climate_SSP370_2090",
pattern = ".tif$",
full.names = TRUE)
# Load all bioclim rasters
future_bioclim <- lapply(files, terra::rast)
# Make into one raster stack
future_bioclim <- rast(future_bioclim)
# Plot the future bioclimatic variables
plot(future_bioclim)class : SpatRaster
dimensions : 681, 841, 19 (nrow, ncol, nlyr)
resolution : 0.05, 0.05 (x, y)
extent : 111.975, 154.025, -44.025, -9.975 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
sources : bioclim_01.tif
bioclim_02.tif
bioclim_03.tif
... and 16 more sources
names : bioclim_01, bioclim_02, bioclim_03, bioclim_04, bioclim_05, bioclim_06, ...
[1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n MEMBER[\"World Geodetic System 1984 (G2296)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"
Code
class : SpatRaster
dimensions : 834, 898, 19 (nrow, ncol, nlyr)
resolution : 5590.925, 5590.925 (x, y)
extent : -4407995, 612655.7, 4234346, 8897178 (xmin, xmax, ymin, ymax)
coord. ref. : GDA2020 / MGA zone 56 (EPSG:7856)
source(s) : memory
names : bioclim_01, bioclim_02, bioclim_03, bioclim_04, bioclim_05, bioclim_06, ...
min values : 0.00000, 0.00000, 0.0000, 0.0000, 0.00000, -2.243844, ...
max values : 33.02249, 16.09955, 66.1461, 706.1523, 46.24717, 25.827158, ...
Mask to SEQ extent
We can see that for these layers the water surrounding Australia is not NA, but are values of 0. In some cells this introduces an artifact on the coastline where the values in some cells are an average between realistic values and 0, which results in cells with unusual and unrealistic values. We will do our best to mask these out.
Code
We can clearly see the artifacts in BIO5, which are not present in the current layers.
If we plot a histogram of the values we can see where we can set a threshold to exclude those lower values on the coastline. We tested values to ensure that we were removing as many of the dodgy values at the coast whilst retaining all real values, and the best balance was 28 degrees.
Unfortunately now we have slightly fewer cells on the coastline, but as these were artifacts we don’t have the true data for those cells anyway.
Code
Code
Code
We can see there are a couple raster cells that have artifacts, but we can filter by the difference between the current and future layers to remove these. Unfortunately, temperature is only going to increase, so we can filter by values where the difference is negative for the BIOCLIM 5 layer, which is the max temperature of the warmest month.
Code
[1] -1.306969
Code
Code
Now that we have applied the artifact mask across all layers, they look better.
Now we can save the future environmental covariates for SEQ.
Load future environmental data 2
Here we load outputs from a low emissions shared socio-economic path scenario (SSP 1.26) for the year 2090 (2080 - 2099).
Code
class : SpatRaster
dimensions : 681, 841, 19 (nrow, ncol, nlyr)
resolution : 0.05, 0.05 (x, y)
extent : 111.975, 154.025, -44.025, -9.975 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
sources : bioclim_01.tif
bioclim_02.tif
bioclim_03.tif
... and 16 more sources
names : bioclim_01, bioclim_02, bioclim_03, bioclim_04, bioclim_05, bioclim_06, ...
min values : 0.00000, 0.00000, 0.0000, ? , ? , ? , ...
max values : 30.47033, 16.75114, 63.2751, ? , ? , ? , ...
[1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n MEMBER[\"World Geodetic System 1984 (G2296)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"
Code
class : SpatRaster
dimensions : 834, 898, 19 (nrow, ncol, nlyr)
resolution : 5590.925, 5590.925 (x, y)
extent : -4407995, 612655.7, 4234346, 8897178 (xmin, xmax, ymin, ymax)
coord. ref. : GDA2020 / MGA zone 56 (EPSG:7856)
source(s) : memory
names : bioclim_01, bioclim_02, bioclim_03, bioclim_04, bioclim_05, bioclim_06, ...
min values : 0.00000, 0.00000, 0.00000, 0.0000, 0.00000, -4.081378, ...
max values : 30.78085, 16.89395, 63.34254, 687.5474, 43.83962, 23.968334, ...
Mask to SEQ extent
We will crop and mask the future bioclimatic variables to the SEQ extent, and apply the artifact mask as we did for the SSP 3.70 layers.
Code
# First mask and crop by the SEQ extent
future_bioclim <- terra::mask(future_bioclim, SEQ_extent.vect)
future_bioclim <- terra::crop(future_bioclim, SEQ_extent.vect)
# set all artifact values to NA
future_bioclim <- terra::mask(future_bioclim, artifact_mask, maskvalues = 0) # Set all values of 0 to NA
future_bioclim <- terra::mask(future_bioclim, temp_diff_mask, maskvalues = 0) # also use the temp difference mask
plot(future_bioclim[[5]])